library(rstan)

cat("\nCCES 2016 PID\n\n")

# Set to the local directory of the replication data
setwd("~/Downloads/Replication_Data/")
n_cores <- 8
options(mc.cores = n_cores)

L2016D <- readRDS("Data/CCES-L2016D.rds")
L2016R <- readRDS("Data/CCES-L2016R.rds")

##### Stack data to feed to Stan

model_data_2016D <- list(n_items = length(unique(L2016D$j)),
                        n_respondents = length(unique(L2016D$k)),
                        n_obs = nrow(L2016D),
                        k = L2016D$k, j = L2016D$j, y = L2016D$y)

model_data_2016R <- list(n_items = length(unique(L2016R$j)),
                        n_respondents = length(unique(L2016R$k)),
                        n_obs = nrow(L2016R),
                        k = L2016R$k, j = L2016R$j, y = L2016R$y)

##### Compile to C

model_am <- stan_model(file = "AM_Sigma.stan")

##### Sample from posterior

exclude_pars <- c("beta_raw", "alpha_raw", "zeta_raw")

# Democrats
n_items <- length(unique(L2016D$j))
inits <- replicate(n_cores, list(list(zeta_raw = c(-1, 1, -1, 1, -1, runif(n_items-7, -2, 2)))))

posterior_cces_2016D <- sampling(model_am, data = model_data_2016D,
                                warmup = 1000, iter = 1250, seed = 1,
                                include = FALSE, pars = exclude_pars,
                                init = inits,
                                refresh = 50, save_warmup = FALSE,
                                control = list(adapt_delta = 0.8,
                                               max_treedepth = 10,
                                               stepsize = 1),
                                open_progress = TRUE, thin = 1, chains = n_cores)

saveRDS(posterior_cces_2016D, "Fitted_Models/posterior_cces_2016D.rds", compress = "xz")


# Republicans
n_items <- length(unique(L2016R$j))
inits <- replicate(n_cores, list(list(zeta_raw = c(-1, 1, -1, 1, -1, runif(n_items-7, -2, 2)))))

posterior_cces_2016R <- sampling(model_am, data = model_data_2016R,
                                warmup = 1000, iter = 1250, seed = 1,
                                include = FALSE, pars = exclude_pars,
                                init = inits,
                                refresh = 50, save_warmup = FALSE,
                                control = list(adapt_delta = 0.8,
                                               max_treedepth = 10,
                                               stepsize = 1),
                                open_progress = TRUE, thin = 1, chains = n_cores)

saveRDS(posterior_cces_2016R, "Fitted_Models/posterior_cces_2016R.rds", compress = "xz")

